Lake Shore Drive District 6 Crash Analysis

Data Loading and Processing

First, we’ll load our dependencies and process the crash data:

# Source utility functions (which loads libraries)
source("R/utils.R")

# Read and process crash data
crashes <- load_crash_data("data/raw/crashes/traffic_crashes.csv")
crashes_sf <- convert_to_sf(crashes)

# Read and process fatality data
crash_fatalities <- load_fatality_data("data/raw/crashes/crash_fatalities.csv")
crash_fatalities_sf <- convert_to_sf(crash_fatalities, "Longitude", "Latitude")

# Combine crashes and fatalities
combined_crashes_sf <- bind_rows(crashes_sf, crash_fatalities_sf)

District and Street Processing

Next, we’ll process the senate district data and Lake Shore Drive information:

# Load senate district shapefile and filter for lakefront districts
senate_districts <- load_district_shapefile(
  "data/raw/shapefiles/senate/senate_districts.shp",
  new_name = "senate_district"
)

Reading layer `senate_districts' from data source 
  `/Users/mmclean/dev/r/ilga_district_crashes/data/raw/shapefiles/senate/senate_districts.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 59 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -91.51308 ymin: 36.9703 xmax: -87.4952 ymax: 42.50848
Geodetic CRS:  NAD83

lakefront_districts <- senate_districts %>%
  filter(senate_district %in% c("6", "7", "13"))

district_6 <- senate_districts %>%
  filter(senate_district == "6")

# Spatial join crashes with lakefront districts
lakefront_crashes <- st_join(combined_crashes_sf, lakefront_districts, left = FALSE)

# Transform to meter-based CRS for accurate distance calculations
lakefront_crashes <- st_transform(lakefront_crashes, 3857)

# Keep district 6 crashes for district-specific analysis
crashes_with_district <- lakefront_crashes %>%
  filter(senate_district == "6")

# Query OSM for Lake Shore Drive
lake_shore_streets <- opq(bbox = CHICAGO_BBOX) %>%
  add_osm_feature(key = 'highway') %>%
  add_osm_feature(key = 'name', value = 'Lake Shore', value_exact = FALSE) %>%
  osmdata_sf() %>%
  pluck("osm_lines") %>%
  {rename(., street_name = name)} %>%
  st_transform(crs = 3857)

Crash Analysis

Process crashes within the Lake Shore Drive buffer:

# Create buffer around Lake Shore Drive and identify intersecting crashes
buffer <- st_buffer(lake_shore_streets, dist = BUFFER_DISTANCE_M)
intersects <- st_intersects(crashes_with_district, buffer, sparse = FALSE)

# Filter crashes within buffer and calculate counts for all lakefront districts
lakefront_intersects <- st_intersects(lakefront_crashes, buffer, sparse = FALSE)
lakefront_crashes$is_lakeshore_crash <- apply(lakefront_intersects, 1, any)
lakeshore_crashes <- lakefront_crashes %>%
  filter(is_lakeshore_crash) %>%
  mutate(
    crash_count = 1,
    fatality_count = ifelse(!is.na(FATALITY_VICTIM) & FATALITY_VICTIM != "", 1, 0),
    total_injuries = INJURIES_TOTAL,
    incap_injuries = INJURIES_INCAPACITATING
  )

# Create fatality subset
lakeshore_fatalities <- lakeshore_crashes %>%
  filter(fatality_count > 0)

# Create district-specific subsets for the time series visualizations
district6_crashes <- lakeshore_crashes %>%
  filter(senate_district == "6")

Crash Impact Visualizations

Fatalities Over Time

# Prepare fatalities data
fatalities_by_year <- district6_crashes %>%
  filter(fatality_count > 0) %>%
  st_drop_geometry() %>%
  mutate(year = lubridate::year(CRASH_DATE)) %>%
  group_by(year) %>%
  summarize(total_fatalities = sum(fatality_count, na.rm = TRUE))

# Get min and max years for consistent x-axis
min_year <- min(fatalities_by_year$year)
max_year <- max(fatalities_by_year$year)

# Create the fatalities visualization
ggplot(fatalities_by_year, aes(x = year, y = total_fatalities)) +
  geom_col(fill = "#9370DB", width = 0.7) +
  geom_text(
    aes(label = total_fatalities),
    vjust = -0.5,
    size = 3.5,
    family = "Arial"
  ) +
  labs(
    x = "Year",
    y = NULL,
    title = "Traffic Fatalities Along Lake Shore Drive in District 6",
    subtitle = "Annual count of traffic-related fatalities",
    caption = "Data: Chicago City Data Portal \"Traffic Crashes - Vision Zero Chicago Traffic Fatalities\" Retrieved Feb. 6, 2025
Includes crashes reported within 100ft of DuSable Lake Shore Drive"
  ) +
  scale_x_continuous(
    breaks = seq(min_year, max_year, by = 1)
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.15)),
    breaks = function(x) seq(0, ceiling(max(x)), by = 1)
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(face = "bold", size = 12, margin = margin(b = 10)),
    plot.subtitle = element_text(size = 10, color = "#666666"),
    plot.background = element_rect(fill = "#FFFCF7"),
    panel.background = element_rect(fill = "#FFFCF7"),
    panel.grid.major.y = element_line(color = "#f0f0f0", linewidth = 0.3),
    text = element_text(family = "Arial"),
    axis.line = element_line(color = "#666666"),
    axis.ticks = element_line(color = "#666666"),
    axis.title = element_text(color = "#666666"),
    axis.text = element_text(color = "#666666"),
    plot.caption = element_text(
      color = "#666666",
      size = 8,
      margin = margin(t = 10),
      hjust = 0
    ),
    plot.margin = margin(t = 15, r = 20, b = 15, l = 15, unit = "pt")
  )

Total Injuries Over Time

# Prepare injuries data
injuries_by_year <- district6_crashes %>%
  st_drop_geometry() %>%
  mutate(
    year = lubridate::year(CRASH_DATE),
    total_injuries = INJURIES_TOTAL
  ) %>%
  group_by(year) %>%
  summarize(total_injuries = sum(total_injuries, na.rm = TRUE))

# Create the injuries visualization
ggplot(injuries_by_year, aes(x = year, y = total_injuries)) +
  geom_col(fill = "#9370DB", width = 0.7) +
  geom_text(
    aes(label = total_injuries),
    vjust = -0.5,
    size = 3.5,
    family = "Arial"
  ) +
  labs(
    x = "Year",
    y = NULL,
    title = "Injuries From Crashes Along DuSable Lake Shore Drive in District 6",
    subtitle = "Total persons sustaining fatal, incapacitating, non-incapacitating, and possible injuries",
    caption = "Data: Chicago City Data Portal \"Traffic Crashes\" Retrieved Feb. 6, 2025
Includes crashes reported within 100ft of DuSable Lake Shore Drive"
  ) +
  scale_x_continuous(
    breaks = seq(min_year, max_year, by = 1)
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.15)),
    breaks = function(x) seq(0, ceiling(max(x)), by = 50)
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(face = "bold", size = 12, margin = margin(b = 10)),
    plot.subtitle = element_text(size = 10, color = "#666666"),
    plot.background = element_rect(fill = "#FFFCF7"),
    panel.background = element_rect(fill = "#FFFCF7"),
    panel.grid.major.y = element_line(color = "#f0f0f0", linewidth = 0.3),
    text = element_text(family = "Arial"),
    axis.line = element_line(color = "#666666"),
    axis.ticks = element_line(color = "#666666"),
    axis.title = element_text(color = "#666666"),
    axis.text = element_text(color = "#666666"),
    plot.caption = element_text(
      color = "#666666",
      size = 8,
      margin = margin(t = 10),
      hjust = 0
    ),
    plot.margin = margin(t = 15, r = 20, b = 15, l = 15, unit = "pt")
  )

Incapacitating Injuries Over Time

# Prepare incapacitating injuries data
incap_injuries_by_year <- district6_crashes %>%
  st_drop_geometry() %>%
  mutate(
    year = lubridate::year(CRASH_DATE),
    incap_injuries = INJURIES_INCAPACITATING
  ) %>%
  group_by(year) %>%
  summarize(total_incap_injuries = sum(incap_injuries, na.rm = TRUE))

# Create the incapacitating injuries visualization
ggplot(incap_injuries_by_year, aes(x = year, y = total_incap_injuries)) +
  geom_col(fill = "#9370DB", width = 0.7) +
  geom_text(
    aes(label = total_incap_injuries),
    vjust = -0.5,
    size = 3.5,
    family = "Arial"
  ) +
  labs(
    x = "Year",
    y = NULL,
    title = "Incapacitating Injuries Along DuSable Lake Shore Drive in District 6",
    subtitle = "Total persons sustaining incapacitating/serious injuries in a crash",
    caption = "Data: Chicago City Data Portal \"Traffic Crashes\" Retrieved Feb. 6, 2025
Includes crashes reported within 100ft of DuSable Lake Shore Drive"
  ) +
  scale_x_continuous(
    breaks = seq(min_year, max_year, by = 1)
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.15)),
    breaks = function(x) seq(0, ceiling(max(x)), by = 5)
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(face = "bold", size = 12, margin = margin(b = 10)),
    plot.subtitle = element_text(size = 10, color = "#666666"),
    plot.background = element_rect(fill = "#FFFCF7"),
    panel.background = element_rect(fill = "#FFFCF7"),
    panel.grid.major.y = element_line(color = "#f0f0f0", linewidth = 0.3),
    text = element_text(family = "Arial"),
    axis.line = element_line(color = "#666666"),
    axis.ticks = element_line(color = "#666666"),
    axis.title = element_text(color = "#666666"),
    axis.text = element_text(color = "#666666"),
    plot.caption = element_text(
      color = "#666666",
      size = 8,
      margin = margin(t = 10),
      hjust = 0
    ),
    plot.margin = margin(t = 15, r = 20, b = 15, l = 15, unit = "pt")
  )

Summary Statistics

# Calculate Lake Shore Drive statistics
lakeshore_total_crashes <- nrow(district6_crashes)
lakeshore_total_injuries <- sum(district6_crashes$INJURIES_TOTAL, na.rm = TRUE)
lakeshore_total_incap_injuries <- sum(district6_crashes$INJURIES_INCAPACITATING, na.rm = TRUE)
lakeshore_total_fatalities <- sum(district6_crashes$fatality_count, na.rm = TRUE)

# Calculate District 6 total statistics
district_total_crashes <- nrow(crashes_with_district)
district_total_injuries <- sum(crashes_with_district$INJURIES_TOTAL, na.rm = TRUE)
district_total_incap_injuries <- sum(crashes_with_district$INJURIES_INCAPACITATING, na.rm = TRUE)
district_total_fatalities <- sum(!is.na(crashes_with_district$FATALITY_VICTIM) & crashes_with_district$FATALITY_VICTIM != "", na.rm = TRUE)

# Create summary statistics dataframe
crash_stats <- tibble(
  Metric = c("Total Crashes", "Total Injuries", "Incapacitating Injuries", "Fatalities"),
  `Lake Shore Drive` = c(lakeshore_total_crashes, lakeshore_total_injuries, lakeshore_total_incap_injuries, lakeshore_total_fatalities),
  `District 6 Total` = c(district_total_crashes, district_total_injuries, district_total_incap_injuries, district_total_fatalities)
) %>%
  mutate(
    `Percentage of District Total` = sprintf("%.1f%%", (`Lake Shore Drive` / `District 6 Total`) * 100)
  )

# Create formatted table
knitr::kable(crash_stats, format = "html") %>%
  kableExtra::kable_styling(
    bootstrap_options = "striped",
    full_width = FALSE,
    position = "left"
  ) %>%
  kableExtra::row_spec(0, bold = TRUE, color = "white", background = "#4B0082") %>%
  kableExtra::column_spec(1, bold = TRUE)

Metric	Lake Shore Drive	District 6 Total	Percentage of District Total
Total Crashes	6248	32490	19.2%
Total Injuries	1415	5452	26.0%
Incapacitating Injuries	166	584	28.4%
Fatalities	16	26	61.5%

Interactive Map Visualization

District 6 Fatality Map

Create an interactive map showing fatality locations in District 6:

# Transform geometries to WGS84 for mapping
lakeshore_fatalities <- st_transform(lakeshore_fatalities, 4326)
lake_shore_streets <- st_transform(lake_shore_streets, 4326)
buffer <- st_transform(buffer, 4326)
district_6 <- st_transform(district_6, 4326)

# Create Lake Shore Drive district 6 fatality map
district6_fatality_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  # Add district 6
  addPolygons(
    data = district_6,
    fillColor = "#9370DB",  # Medium purple
    fillOpacity = 0.3,
    color = "#4B0082",  # Indigo
    weight = 2,
    label = ~paste("Senate District", senate_district)
  ) %>%
  # Add Lake Shore Drive
  addPolylines(
    data = lake_shore_streets,
    color = "#CCCCCC",
    weight = 3
  ) %>%
  # Add the buffer area
  addPolygons(
    data = buffer,
    fillColor = "#CCCCCC",  # Light gray
    fillOpacity = 0.2,
    color = "#CCCCCC",  # Light gray
    weight = 2
  ) %>%
  # Add clustered fatality circles
  addCircleMarkers(
    data = lakeshore_fatalities,
    lng = ~st_coordinates(geometry)[, 1],
    lat = ~st_coordinates(geometry)[, 2],
    radius = 6,  # Base size for individual points
    color = "red",
    fillColor = "red",
    fillOpacity = 0.6,
    stroke = TRUE,
    weight = 1,
    popup = ~paste(
      "Date:", CRASH_DATE, "<br>",
      "Fatalities:", fatality_count, "<br>",
      "Type:", FIRST_CRASH_TYPE
    ),
    clusterOptions = markerClusterOptions(
      maxClusterRadius = 50 * 0.0003,  # Approximate conversion of 50ft to degrees
      spiderfyOnMaxZoom = TRUE,
      iconCreateFunction = JS("
        function(cluster) {
          var count = cluster.getChildCount();
          var size = Math.min(50, 20 + Math.sqrt(count) * 8);
          return new L.DivIcon({
            html: '<div style=\"background-color: rgba(255, 0, 0, 0.6); width: ' + size + 'px; height: ' + size + 'px; margin-left: -' + (size/2) + 'px; margin-top: -' + (size/2) + 'px; border-radius: 50%; display: flex; align-items: center; justify-content: center; color: white; font-weight: bold; border: 2px solid rgba(255, 0, 0, 0.8);\">' + count + '</div>',
            className: 'marker-cluster-custom',
            iconSize: L.point(size, size)
          });
        }
      ")
    )
  ) %>%
  setView(lng = -87.6298, lat = 41.8781, zoom = 12)

# Display the map
district6_fatality_map

Lake Shore Drive Injury Map

Create interactive maps showing injuries, incapacitating injuries, and fatalities along Lake Shore Drive:

# Transform geometries to WGS84 for mapping
lakeshore_crashes_4326 <- st_transform(lakeshore_crashes, 4326)
lakefront_districts_4326 <- st_transform(lakefront_districts, 4326)

# Create binary indicators for crashes with injuries
lakeshore_crashes_4326 <- lakeshore_crashes_4326 %>%
  mutate(
    has_injuries = ifelse(INJURIES_TOTAL > 0, 1, 0),
    has_incap_injuries = ifelse(INJURIES_INCAPACITATING > 0, 1, 0)
  )

# Create injuries heatmap with adjusted parameters
injuries_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  # Add lakefront districts outlines
  addPolygons(
    data = lakefront_districts_4326,
    fillOpacity = 0,
    color = "#4B0082",  # Indigo
    weight = 2,
    label = ~paste("Senate District", senate_district)
  ) %>%
  # Add Lake Shore Drive
  addPolylines(
    data = lake_shore_streets,
    color = "#666666",
    weight = 3
  ) %>%
  # Add injuries heatmap
  addHeatmap(
    data = lakeshore_crashes_4326 %>% filter(has_injuries == 1),
    lng = ~st_coordinates(geometry)[, 1],
    lat = ~st_coordinates(geometry)[, 2],
    intensity = ~1,  # Each crash with injuries counts as 1
    blur = 15,
    max = 0.1,
    radius = 12
  ) %>%
  setView(lng = -87.6298, lat = 41.8781, zoom = 12)

# Display the map
injuries_map

Lake Shore Drive Incapacitating Injury Map

# Create incapacitating injuries heatmap
incap_injuries_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  # Add lakefront districts outlines
  addPolygons(
    data = lakefront_districts_4326,
    fillOpacity = 0,
    color = "#4B0082",  # Indigo
    weight = 2,
    label = ~paste("Senate District", senate_district)
  ) %>%
  # Add Lake Shore Drive
  addPolylines(
    data = lake_shore_streets,
    color = "#666666",
    weight = 3
  ) %>%
  # Add incapacitating injuries heatmap
  addHeatmap(
    data = lakeshore_crashes_4326 %>% filter(has_incap_injuries == 1),
    lng = ~st_coordinates(geometry)[, 1],
    lat = ~st_coordinates(geometry)[, 2],
    intensity = ~1,  # Each crash with incapacitating injuries counts as 1
    blur = 15,
    max = 0.1,
    radius = 12
  ) %>%
  setView(lng = -87.6298, lat = 41.8781, zoom = 12)

# Display the map
incap_injuries_map

### Lake Shore Drive Fatality Map

# Create fatalities map with clusters
fatalities_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  # Add lakefront districts outlines
  addPolygons(
    data = lakefront_districts_4326,
    fillOpacity = 0,
    color = "#4B0082",  # Indigo
    weight = 2,
    label = ~paste("Senate District", senate_district)
  ) %>%
  # Add Lake Shore Drive
  addPolylines(
    data = lake_shore_streets,
    color = "#666666",
    weight = 3
  ) %>%
  # Add clustered fatality circles
  addCircleMarkers(
    data = lakeshore_fatalities,
    lng = ~st_coordinates(geometry)[, 1],
    lat = ~st_coordinates(geometry)[, 2],
    radius = 6,
    color = "red",
    fillColor = "red",
    fillOpacity = 0.6,
    stroke = TRUE,
    weight = 1,
    popup = ~paste(
      "Date:", CRASH_DATE, "<br>",
      "Fatalities:", fatality_count, "<br>",
      "Type:", FIRST_CRASH_TYPE
    ),
    clusterOptions = markerClusterOptions(
      maxClusterRadius = 50 * 0.0003,
      spiderfyOnMaxZoom = TRUE,
      iconCreateFunction = JS("
        function(cluster) {
          var count = cluster.getChildCount();
          var size = Math.min(50, 20 + Math.sqrt(count) * 8);
          return new L.DivIcon({
            html: '<div style=\"background-color: rgba(255, 0, 0, 0.6); width: ' + size + 'px; height: ' + size + 'px; margin-left: -' + (size/2) + 'px; margin-top: -' + (size/2) + 'px; border-radius: 50%; display: flex; align-items: center; justify-content: center; color: white; font-weight: bold; border: 2px solid rgba(255, 0, 0, 0.8);\">' + count + '</div>',
            className: 'marker-cluster-custom',
            iconSize: L.point(size, size)
          });
        }
      ")
    )
  ) %>%
  setView(lng = -87.6298, lat = 41.8781, zoom = 12)

# Display the map
fatalities_map